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ABSTRACT 

Numerical tables available for M/E^/c queueing systems 
are discussed. A new approximation method for steady— state 
information and waiting time distribution of this queueing 
system was developed. Validity of approximation was estab- 
lished directly for the large waiting times and by simula- 
tion for the smaller values. The developed method enables 
one to find delay probability, expected number in the 
queue and in the system, expected time to be spent in the 
queue and in the system, and probability of waiting for 
more than a specified time t. 
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I. 



INTRODUCTION AND SUMMARY 



Multi— server queues constitute a large proportion of 
queueing systems which arise in practice. Among those, 
queues with Poisson arrivals and Erlang service times (which 
will be denoted as M/E^/c) occupy an important position, as 
seen in the banks, airport check— in counters, hotels, super- 
markets etc. However, no significant theoretical work was 
done until the late 60* s. In recent years, methods for the 
steady state solution of M/E^/c queues became available, fol- 
lowed by the numerical tables which were obtained by imple- 
mentation of these methods. Even so, there is still some 
computational difficulty involved obtaining some particular 
information about above-mentioned queues, such as probabili- 
ty of delay exceeding a specified time length. 

The importance of the M/E y/c queue is due to the fact 
that it is used to model any queueing system whose service 
time distribution is believed to be unimodal. The solutions 
are known for extreme values of k: exponential service time 
distribution for k=l and constant service times as k-*-°°. 

Usable solutions for multi— server queues having either expo- 
nential or constant service times are readily available. 

The M/E^/c queue is also important by providing a large 
variety of service time distributions in between these two 
extremes . 
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The purpose of this study is, by means of computer simu- 
lation and examination of numerical tables published earlier, 
to present a simple and accurate computation method for ob- 
taining information about steady— state solution of the M/E^/c 
queue. 

For this purpose, a general description of M/E^/c queue- 
ing model with definitions and assumptions for the analysis 
of it are given in Chapter II, followed by a summary of for- 
mer studies. 

In Chapter III, the numerical tables are analyzed and 
some conclusions are drawn in terms of a simple form approxi- 
mation for steady— state distribution of waiting times. 

This approximation for the distribution of’ waiting times 
needed verification for small values of them. These cases 
were studied by computer simulation. The procedure and the 
results of simulation are given in Chapter IV. 

The results of the analysis of numerical tables and simu- 
lation are combined and generalized in Chapter V. Then the 
computational method developed by this generalization is de- 
scribed on a step— by— step basis. With this method, there is 
no need for numerical tables, computations are very simple, 
and the results are accurate for most practical purposes 
(the error being about 2% in general) . To demonstrate the 
advantages of the proposed method, a comparison of it with 
the other method which uses numerical tables is given at the 
end of that chapter. 
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II. The M/E j,/ c Queue 



The description of the M/E-^/c queueing model with assump- 
tions, definition of parameters and system variables along 
with some basic definitions in the queueing theory forms the 
first half of this chapter. A brief discussion of the pre- 
vious studies on this model then follows. 

A. DESCRIPTION OF THE SYSTEM 

The M/Ej,/c queue is a multi— server queueing system with 
c servers, where arrivals occur according to a Poisson pro- 
cess with rate X, and service times have Erlang distribution 
with shape parameter k. It is also an infinite queue, i.e. 
there is no limitation for the number of customers in the 
system. 

The distribution of interarrival times X is given by 



and the service times have the following probability density 
function: 



f X (x) = 



— Xx 



x > 0 



f 



f Y (y) 




/ 



y > 0 



(k— 1) ! $ 



2 

with mean kg and variance k$ 
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1 . Definitions of Parameters 

a. Arrival Rate and Service Rate 

Arrival rate to the system is the reciprocal of 
mean interarrival time and denoted by X. Similarly, service 
rate of a server is defined as reciprocal of mean service 
time and represented by y. Then, y = . 

Kp 

b. Traffic Intensity 

Traffic intensity p is defined by the following 

expression: 



p = X/cy . 

c. Offered Load 

Offered load is the ratio of arrival rate X to 
service rate y and denoted by a. It can also be expressed 
by product of number of servers c and traffic intensity p 
as follows 



a = cp . 



d. Coefficient of Variation 

The ratio of standard deviation to mean is de- 
fined as coefficient of variation of a probability distribu- 
tion. It is denoted by V. For Erlang distribution coeffi- 
cient of variation is expressed as 



V 



/ke 2 

kB 



_ 1 _ 

/k 
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2 . System Variables 



a. Number in the System 

The total of the .number of customers waiting in 
the queue and those in the service is defined as number in 
the system and denoted by N. The number of customers only 
in the queue is represented by Nq. Expected values of N 
and Nq are denoted by L (average number in the system) and 
Lq (average number in the queue) , respectively. 

b. Waiting Time 

The time spent in the queue by a customer is de- 
fined to be the waiting time of a customer. The time spent 
in the system is then the sum of waiting time and service 
time of a customer. Tq denotes waiting time and T denotes 
the time spent in the system. Expected values of T and Tq 
are denoted by W (average time in the system) and Wq (aver- 
age waiting time) respectively. 

c. State Probabilities 

The number in the system at a particular time 
is defined as the state of the system. Then state probabili- 
ty P (N = n|t) is the probability of the system being in 
state n at time t. 

d. Delay Probability 

Delay probability is the probability that the 
arriving customer finds all the servers busy and enters the 
waiting line. It is denoted by P (Tq >0). Also, by defini- 
tion it follows that 



P (Tq > 0) = P(N _> c) . 
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3 . The Steady— State Solution 



The steady— state solution is defined to be the proba- 
bility distribution of the number in the system when the 
system achieves statistical equilibrium. In the steady— state, 
the state probabilities do not depend on t, i.e. P(N = n|t) = 
P (N = n) . Also it follows that 



X P(N = n) = 1 . 
n=0 

In queueing theory, it is a well-known fact that [1] 
the steady— state is achieved if and only if the traffic inten- 
sity p is less than 1. 

4 . Assumptions 

The M/E^/c queueing system under study is assumed to 
have the following properties: 

(i) There is only one waiting line no matter how 
many servers there are. 

(ii) The customer at the head of the waiting line 
starts getting service immediately whenever a server becomes 
available. 

(iii) Order of service is first— come first— served. 

(iv) If an arriving customer finds all servers busy, 
he joins the waiting line. The waiting line has unlimited 
capacity. 

(v) The service channels (servers) are indexed 
consecutively. If an arriving customer finds more than one 
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server vacant, he goes to the Server with smallest index 
(this assumption doesn't cause any loss of generality, but 
is convenient for the computer simulation model) . 

(vi) The servers are homogeneous, i.e. the service 
distribution is the same for all servers. 

(vii) Only one customer at a time can be served by 
a server. 



B. FORMER STUDIES 



There are numerous studies in the literature of queueing 
theory, done for the steady— state solution of M/E^/c queue. 
But most of them cover only some particular values of k and 
c, and their results cannot be generalized. Therefore, such 
studies will not be mentioned here individually. 

The earliest suggestion known by the author about the 
solution of general M/E ,/c system was mentioned by Lee [2] . 
He stated referring to a case study done in 1956 that mean 
waiting time can be approximated by use of David Owen's sug- 
gestion. The formula given in [2] was said to be usable for 
0.7 < V < 1.0 and as follows. 



W q - I W ql (1+V2 > 



where 



W 



’ : mean waiting time for M/E^/c queue 

W mean waiting time for M/M/c queue 

V : coefficient of variation of Erlang (k) 
distribution, = 1/k. 

The validity of the approximation was demonstrated by 
comparing its results with simulation results for different 
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cases. However, the approximation is just for mean waiting 
time or queue length, hence it isn't possible to approximate 
state probabilities. Also, no reference was given for de- 
tails of David Owen's suggestion. This approximation will 
be mentioned again later in this study. 

The steady— state solution of the general M/Ej,/c queue 
was first studied by Mayhugh and McCormick [3] . The results 
can be applied to the cases with any value of k and c. How- 
ever, the computation procedure is so complex that it 
requires a very considerable amount of work. 

A parallel study was also done by Heffer [4] . The solu- 
tion method proposed by Heffer differs from Mayhugh and 
McCormick's, but it also is very complex. 

The reader is referred to the references [5] and [6] for 
more detailed discussion of these two studies. 

The most recent study was done by Yu [5] . It concerned 
the E m /E^/c queue with heterogeneous servers, but results 
were also stated for homogeneous server case. Then letting 
m = 1 gives the solution procedure for M/E^/c queue. However, 
like the other two studies, the computations required still 
demand an enormous amount of work. 

The theoretical results obtained by Heffer [4] and Yu [5] 
formed the basis of the only numerical tables available for 
M/E^/c queue, prepared by Hillier and Lo [6] . The tables 
have cases of M/E^/c queue with limited values of k and c 
(k < 0, c < 10), and a few cases for E /E v /c queue. These 
tables will be discussed in detail in the following chapter. 
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III. DISCUSSION OF NUMERICAL TABLES AND 
DEVELOPING A NEW COMPUTATION METHOD 



In this chapter, a detailed review and discussion of the 
numerical tables given in [6] will be made. Then, using the 
results of these discussions, a new computation method will 
be developed. 

Sample pages from the numerical tables are given in Appendix B. 

A. DISCUSSION OF NUMERICAL TABLES 
1 . General Description 

The tables under consideration were designed for gen- 
eral E m /E^/c systems. The analysis however will be focused 
on the cases in which m = 1 (i.e. M/E^/c system) , which forms 
a large porportion of the tables. As mentioned earlier,, the 
values of k and c for various tables are as follows: 



k=2 ; c=l,2 



10 



k=3 ; c=l,2 



• • • • f 



5 



k=4 ; c=2 , 3 



k= 5 , 6 , . . . , 8 ; c= 2 



The information given in each table is: 

(i) State probabilities, P(N=n), and 

(ii) Cumulative state probabilities, P(N<n), for 



n=l , 2 



t • • • / 



(iii) Delay probability, P(Tg>0); 

(iv) Expected number in the system, L; 
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(v) Average queue length, Lq; 

(vi) Average queue length for the equivalent M/M/c 



system with the same arrival and service rates, L qi* 

(vii) The ratio Lg/L 

All this information is given in each case and for 
the following values of traffic intensity p : 

p = 0.10, 0.20,..., 0.50, 0.55,..., 0.95, 0.98, 0.99. 

So, all the information is given for steady— state con- 
ditions . 

Average queue length for the equivalent M/M/c sys- 

tem is computed by using the following formula [7] , 

c 



(cp) 



ql , , . .2 o 

M c! (1-p) 



(3.1) 



where 



c— 1 



- I 



(cp) 



(cp) 



j=0 



j ! + c! (1-p) 



The rest of the information given in the tables was 
computed by the results of theoretical work mentioned in the 
last chapter, done by Heffer and Yu. 

2 . Use of Tables 

The kind of information listed above can be obtained 
directly from the tables. However, by making use of some 
general relationships in the queueing theory, some other in- 
formation can be obtained besides those mentioned above, 
a. Interpolation 

No interpolation method is specified in the ex- 
planation sections for the tables in [6] , for the values of 
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p different from those given above. The default position 
taken is linear interpolation, 
b. Waiting Times 

The average waiting time W and average time in 

Si 

the system W are frequently of interest in the analysis of 
queueing systems. These values can be obtained from the 
tabulated values of L and L respectively, by using the fol— 
lowing relationships [8] , 




The probability P (T >t) that waiting time exceed— 

91 

ing t is also important in the analysis of queueing systems. 

Although this kind of information is not given in the tables, 

P(T >t) can be approximated from the state probabilities, 

'“d 

P (N=n) , as follows: 

00 

P (T >t) = 2 P(N-n) P{D<_[k(n-c+l)-l] } , t>0, 

^ n=c 

(3.2) 

where D is a Poisson random variable with mean ckyt. The 
reader is referred to the discussion in [6] for the stochas- 
tic reasoning of this relationship. The quantity 
P{D£[k (n— c+1)— 1] } can be obtained from Poisson distribution 
tables with mean ckyt. If this mean exceeds 25, then the 
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normal distribution with the same mean and variance can be 
used as an approximation to Poisson distribution. 

3. Discussion 

a. Delay Probabilities 

Similar to M/E^/c, M/M/c denotes a multi— server 
queueing system with Poisson arrivals of rate A and exponen- 
tial service times of service rate y. The delay probabili- 
ties for M/M/c queue can be computed by using the formula 
[7], 



P (T g >0) 



(cp) C 
c! (1-p) 



( 



c=l 

V 

j=o 



(cp) 



(cp) c \ 
c! (1-p) / 



(3.3) 



A chart is also available in Appendix A for c=2, 
3,... ,15 and 0<p<l, for P(T >0). 

A comparison of delay probabilities of M/Ej,/c 
and M/M/c systems for selected values of p, c and k is given 
in Table I of this study. A careful examination of this 
table shows that the corresponding delay probabilities of 
the two systems for the same value of p are very close. 

Since the formula or charts are available for M/M/c system, 
this comparison shows that they can also be used to estimate 
delay probabilities for M/E^/c system with a very small 
error, usually less than 0.01. 
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Comparison of Delay Probabilities in Different Systems 
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b. Average Queue Length and Number in the System 
A careful examination of the "Ratio" column of 
Table I in [6] indicates that 



lim _ k+1 _ 1,, 1« 

p-K L L ql 2k k' 



|(1 + V 2 ) 



(3.4) 



which is essentially the same as Owen's suggestion [2], How- 
ever, it works beyond the limits of coefficient of variation 
V that were stated by Owen. Also (3.4) is true even when p 
is markedly less than one. The theoretical verification for 
(3.4) was developed by Yu in [9]. 

A more precise approximation formula was also 
developed by Hillier and Lo in [6] which provides more accu- 
racy for small values of p. In this case, the approximation 
formula is given by 



L 

q 



i(l + g) (1 ♦ £> L ql 



(3.5) 



where g = f(p,k,c). Then the expression for g has been ob- 
tained by linear regression. The exact coefficients for 
g(p,k,c) are given in [6]. However, a rougher expression 
which is more convenient for computations will be used here, 
given as 



g = ( c_1 ) 2/3 Ul-P) + (1-p) 2 } (3.6) 
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The main purpose of Hillier and Lo for introduc- 
ing g into approximation formula (3.5) was to provide means 
of approximation for the cases with larger k and c values 
which were not covered in [6] , namely, extrapolation for lar- 
ger values of c and k. One way to check the validity of the 
computational formula for g is to first let k go to infinity, 
then compare the values of average waiting time computed 
by the formula given below with average waiting time W^ D of 
M/D/c queue. This latter queue has Poisson arrivals with 
rate X and constant service times of length 1/y (or kg) , and 
one can use the fact that the distribution of service times 
can be approximated by constant service time distribution 
for very large k as mentioned in Chapter I. 

The computation formula for W^ D of M/D/c queue 
is given in [2] as follows. 




CO 



I 



— ia 
e 



T (ia) j 

-i ! 



1 

P 



5" (ia) ^ 

^ i 1 

j=ic+l * 



However, a more convenient form for computations can be ob- 
tained as 



W 



qD 



- I 

i=l 



—ia . . >ic 
e (ia) 



(ic) ! 



+ (1 - 



F> S 

p j=ic+l 



e la (ict) 



(3.7) 
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which permits the use of Poisson distribution tables. But 
it still requies a great deal of computational work. For- 
tunately, a very convenient chart was developed by Shelton 
[10] for Wg D for l<c<100 ,and 0.10£p<0. 96. 

The comparison mentioned above can be made hav- 
ing the necessary tools available and Table II can be used 
for this purpose. The values of W^ D for M/D/c system were 
obtained from Shelton's charts. Average waiting time for 
M/E^c system is computed by using the following formula. 



w q - Tit 1 + Tj'E+r’ <c-l> 2/3 Ul-P) + <l-p) 2 >]u + |)L ql 

where L ^ is computed from (3.1) . If the expression given 
by (3.6) is true for very large values of k, then it should 
also be true that 



lim W q = W qD 
k->°° 



Taking the limit of gives 



lim W g = + yj (c-l) 2//3 {(1-p) + (l-p) 2 }J L 



ql 



k-*-°° 



On the other hand, if g is omitted, the limit 

will then be 

lim W = i L . = 4- W . 
q 2X ql 2 ql 

k-*» 
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A comparison 


of W . and these two limits for 
qD 


eight different cases is 


given in 


Table II. 


Investigation 








TABLE 


II. 




Comparison of 


v and 


1 Limits 


of W With 

q 


or Without g 


Case No 


c 


p 


1 W . 

2 ql 


W (k=«>) 

q 


W _ 
qD 


1 


2 


0.80 


1.185 


1.209 


1.08 


2 


2 


0.80 


1.422 


1.451 


1.296 


3 


5 


0.90 


2.287 


2.34 


2.16 


4 


5 


0.90 


1.144 


1.17 


1.08 


5 


10 


0.90 


1.337 


1.391 


1.36 


6 


10 


0.90 


1.505 


1.563 


1.53 


7 


20 


0.90 


3.305 


3.766 


3.36 


8 


20 


0.90 


5.508 


5.88 


5.6 


of these 


results 


shows that, for 


almost all 


the cases the 



limit with g omitted is closer to W . than the other limit. 

qD 

This indicates that the formula (3.6) should be reviewed 
for large k values. However, this revision can be done only 
after some additional tables become available for the cases 
not covered in [6] . 

c. Waiting Times 

The method of approximating the probability 
P (T >t) that the waiting time will exceed t using the state 
probabilities P (N=n) was described earlier in the chapter. 
However, this method assumes that an arriving customer will 
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have to wait for approximately k(n— c+1) service phase (ex- 
ponentially distributed with mean $) completions before a 
server becomes available to start his service. It is stated 
in [6] that this approximation for P (T >t) should be better 
for larger values of p and t due to this assumption. It 
should be added that, for large k and small t, the approxi- 
mation fails greatly in representing the actual distribution 
of waiting times for the same reason. 

B. NEW COMPUTATION METHOD 

The construction of a practical computation method using 
the facts obtained from the preceding discussions will form 
the rest of the chapter. This new method is desired to be 
applicable to the queueing problems involving M/E^/c systems 
met in practice, with simple and quick computational work 
without any need for tables, interpolations etc. 

The kinds of information which are desired to be able to 
compute by the new method for M/E^/c queueing system are 
ma inly 

(i) Delay probability, 

(ii) Average waiting time, average time in the system, 
average queue length and average number in the system, 

(iii) Probability distribution of delay times and proba- 
bility that waiting time will exceed a specified time length. 

1. Delay Probability 

In the discussion earlier it was pointed out that 

the delay probabilities P(T >0) for the systems M/E,/c and 

4 
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M/M/c are very close for the same value of p. An examina- 
tion of the comparisons given in Table I shows that the 
delay probabilities are very close for M/M/c (k=l) , M/E^/c 
(k>l) and M/D/c (k=°°) systems with the same value of p. 

The differences between the delay probabilities for the same 
p are less than 0.01 for almost all the cases and gets even 
smaller for p close to one. This comparison suggests that 
delay probabilities of the M/M/c system can be used for cor- 
responding M/Ej,/c system also. These delay probabilities 
are computed by using formula (3.3). The formula is not 
suitable for large values of c, however, and a chart is given 
in Appendix A for up to about 16. For greater values of c, 
the charts are also available in [10] and [11] . 

2 . Average Queue Length and Average Waiting Time 

In the previous discussion, after examination of the 
tables in [6] it was concluded that the average queue length 
for M/Ej,/c system can be approximated by 

L q = i(l + g) (1 + |)L ql (3.8) 



where L ^ is computed according to (3.1), and g is computed 



by 



V? V 



2 k— 1 , . ; 2/3 , ... , 2-, 

g = i2 k+T (c-1) {(i-p) + (1-p) } 



Contribution of g in equation (3.8) is negligible 
for practical purposes in most cases. Also it was shown 
earlier in this chapter that the formula g requires some 
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I 



modifications. Therefore it will be omitted in further 
discussions . 

Now, attention will be focused on equations (3.1) 
and (3.3). Rewriting equation (3.1) in a different form, 
then substituting (3.3) gives 



ql (1-P) 



(cp) 



c! (1-p) 



rc-i 

£ 

j=o 



(cp) j (cp) c 

3 ! c! (1-p) 



and 



p*P (T g >0) 
ql = (1-p) 



Then it follows that 



L . P (T >0) 

W = _2i = <3 

ql X cp (1-p) 



Using (3.8) with g=0, some approximation formulas 



for average queue length and average waiting time are ob- 
tained respectively as follows 



L_ = 



W = 



1 P(T >0) 

2 p (1-p) 



1 

1 



P (T >0) 
SI 

cy (1-p) 



a ♦ |) 



a + £> 



(3.9) 

(3.10) 



where P(Tg>0) is delay probability for M/M/c system which 
can be obtained from the charts very easily. Then it is 
very simple to compute L and W using equations (3.9) and 

q q 

(3.10) . 
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The average number in the system L and average time 
in the system W can be computed by following well-known rela- 
tionships of queueing theory. 



L = 





3 . Waiting Time Distribution 

Suppose that the conditional waiting time distribu- 
tion given T >0 for M/E, / c queueing system can be approxi— 
q k 

mated by 



P(T g >t 



Tq>0> - 



— bt 



Then using Bayes' theorem it follows that 



P(T >t) = P (T >t | T >0) P (T >0) = e bt P (T >0) . (3.11) 



Let be the CDF of waiting times. Then 



F T q ( t ) = l-P(T g >0) e bt , 



t 0, 



Now average waiting time can be computed as follows 

00 00 

W q = E(T g ) = o f (l-F Tq (t))dt = o / P(T q >0)e” bt dt 



P (T >0) 

w = 2 

q b 



(3.12) 



Suppose the parameter b is modeled as 

t. _ 2cu (1— p) 

D 2 
(1+V Z ) 
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Then , 

P (T >0) - 

W q - Scufl-p) < 1+V > (3 - 13) 

where V is coefficient of variation of service time distri- 
bution. Notice that equation (3.13) is exactly the same as 
equation (3.10). However, there can exist some other wait- 
ing time distribution to give the same average waiting time 
as given by (3.13), so it must be shown that for M/E^/c sys- 
tem, the distribution of waiting times can be approximated 
by 



- 2cy (1-p) t 

F (t) = l-P(T g >0)e ( 1+V 2 ) , t>0 (3.14) 

If (3.14) is true, then P (T >t) computed by (3.11) 
should be approximately equal to the one obtained by the pro- 
cedure suggested by Hillier and Lo [6] which uses state pro- 
babilities as described earlier by equation (3.2). Table 
III gives the comparison of P(T^>t) values obtained by these 
two different formulas for various values of t. It should 
be recalled that approximation by state probabilities is 
poor for small values of t. However, for t>W the two re- 
sults agree very well. The difference is ease of computa- 
tions. Equation (3.11) is much easier than equation (3.2) . 
to compute the same value. 

The validity of approximation (3.14) for small values 
of t will be shown by simulation, which forms the next chap- 
ter. Then the suggested computation method will formally 
be described on a step— by— step basis in last chapter. 
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Comparison of the Two Approximation Methods 
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IV. COMPUTER SIMULATION 



A computer simulation model for M/E^/c system is used 
to investigate the validity of approximation for waiting 
time distribution introduced in the preceding chapter. 
Validity of this approximation for moderate and large 
values of t was demonstrated in the same chapter. There- 
fore, analysis will now be focused on the approximation 
for small values of t. 

A. COMPUTER PROGRAM 
1. Model 

Since the model used to construct the simulation 
program is based on the same assumptions as stated in Chap- 
ter II, it won't be described again in this chapter. The 
computer language selected was SIMSCRIPT which is especial- 
ly convenient for simulation of multi— server queueing sys- 
tems. The reader is referred to [12] and [13] for a brief 
idea about SIMSCRIPT if he is not familiar with this 
language. Main reference for SIMSCRIPT however is [14] . 

One way to check for small t the validity of approxi- 
mation distribution developed earlier is to keep the fre- 
quency table of waiting times of the customers during the 
simulation. Let M(t) be the number of customers with wait- 
ing time greater than t and M be total number of customers 
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served during simulation period. Then the estimator for 
probability of waiting time greater than t is 



This value can be compared with the value given by the ap- 
proximation formula to check its validity. 



written in the way that it would be possible to obtain the 
estimates of all kinds of information about the queueing 
system being simulated; average waiting time and average 
queue length, state probabilities, delay probability etc. 
However, only the estimates given by (4.1) will be discussed 
here. 

2 . Variance Reduction 

One of the problems encountered in the simulation 
of queueing systems is the high variability of waiting times, 
especially when traffic intensity p is high. Also strong 
positive correlation between the waiting times of consecu- 
tive customers causes serious errors if the standard for- 
mula is used to estimate the sample variance of waiting 
times since this formula will underestimate the true vari- 
ance [1] . 



difficulty with high variability. Interested reader can find 
comprehensive information in references [12] , [15] and [16] . 

The variance reduction method which will be used here is 
antithetic variates. 




(4.1) 



The computer program is given in Appendix B. It is 



There are several methods developed to overcome this 
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In the simulation program, two different random num- 
ber streams are used to generate interarrival and service 
times. Once a uniform variate between 0 and 1 is generated, 
then the random variate from a desired distribution can be 
obtained by 

X = F — 1 (u) (4.2) 

where u is a uniform variate between 0 and 1 and F ^ denotes 

x 

the inverse of CDF of X [12] . Then, for example, a sequence 
of exponentially distributed random variates with mean 1/A 
can be obtained from a stream of uniform variates by using 

X = — j ln(u) (4.3) 

Let and Z 2 be estimators of a parameter, obtained 
from two different simulation samples and Z^ = (Z^ + Z 2 ) • 

Then 

Var(Z 3 ) =ivar(Z 1 +Z 2 ) = j [Var (Z^ +Var ( Z 2 ) +Cov (Z 1 , Z 2 ) ] 

(4.4) 

which implies that maximum negative correlation between Z^ 
and Z 2 would minimize Var(Z 3 ). This can be obtained by us- 
ing 1— u in (4.2) in place of u for random variate genera- 
tion in the second simulation run. 

During the simulation, a random sequence of large 
service times or short interarrival times would cause long 
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waiting times and conversely. Using the procedure described 
above, a sequence of large waiting times would become sequence 
of short waiting times in the second run? in other words, the 
waiting times in two different runs will be negatively corre- 
lated. Let P^(t) be the probability that waiting times will 
exceed t in the ith run, i = 1,2. Using antithetic variates, 
P^(t) and P 2 (t) w iH be negatively correlated; then the esti- 
mator will be computed by 

P(T >t) = j(P 1 (t) + P 2 (t)) 

« 

where P^ (t) and P 2 (t) are computed according to (4.1). 

Generation of antithetic exponential variates for in— 
terarrival times is given by (4.3). However, it is not 
straightforward to generate antithetic Erlang variates since 
the CDF cannot be inverted in closed form. Generating Erlang 
variates as sum of k exponential variates with mean 3 does 
not help since equation (4.2) cannot be utilized to obtain 
antithetic variates. Nevertheless, one way to solve this 
problem is to store the CDF table of chi— square distribution 
in an array, then compute F ”~^(u) by linear interpolation to 
get a chi— square random variate, and obtain the Erlang 
variate by the following transformation 




where E^ is the Erland variate with mean k0 and C^ is the 
chi-square variate with degrees of freedom k. The tables 
for CDF of chi-square distribution are available in [17] . 
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3. Data Collection 



a. Selection of Input Parameters 

The input values for each pair of runs were arbi- 
trarily selected to give certain traffic intensity p values, 
usually 0.8 and 0.9 in order to get a wider range of waiting 
times even though they cause high variability. The width of 
frequency intervals was chosen to be 0.25 min. so that f^, 
for example, would show the number of customers with waiting 
time less than 0.25 minutes. Then M(t) , the number of 
customers with waiting time greater than t can be computed 
by 

00 

M(t) = 2 f • 
i=4t+l 1 

since t will be in multiples of 0.25. This computation is 
done in the computer program. 

b. Initial and Final Conditions of the System 
Initial and final conditions of the simulated 

system are important since they effect the value of para- 
meter estimated. The approach taken in this study was to 
use epochs. If the time at which the number in the system 
N changes from 0 to 1 by an arrival is defined as a regenera- 
tion point, then the interval between two successive regene- 
ration points can be defined as an epoch [12] . It is 
obvious that the sequences of waiting times in two different 
epochs will be independent from each other. The epochs tend 
to be lengthy with high traffic intensity p, large arrival 
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rate A, with large number of servers c, or a combination of 
these three factors. Before each simulation run, the number 
of epochs for which the simulation is to be run was deter- 
mined as an input considering those three factors. Exper- 
ience showed that the first few epochs tend to have waiting 
times smaller than usual, so they were omitted for data 
analysis. This isn't feasible for cases with long epochs; 
however the first couple of epochs have enough information. 
The simulation runs were started with N = 0 and an arrival 
so that starting time would be a regeneration point and 
stopped after the number of epochs determined earlier is 
completed, i.e. N = 0 was the final condition. 

P^(t) or P£(t) are given in the output of program 
for t = 0.25, 0.50, 0.75, .... etc. Then to get the estima- 
tors P(Tg>t), all one has to do is to average them for cor- 
responding t values. 

B . RESULTS 

A comparison of waiting time probabilities obtained from 
simulation and approximation method is given in Table IV. 
Investigation of this table indicates that the approximation 
method agrees quite well with the results of simulation for 
small t values. The difference between the corresponding 
values for moderate or large values of t for which the 
validity of approximation was demonstrated earlier explains 
the difference between the values when t is small. 
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Comparison of Simulation and Approximation Results 
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V. CONCLUSION 



\ 



It was shown in the last chapter that the approximation 
method developed earlier for the distribution of waiting 
times can also be used for small values of t. The method 
for computing steady— state information of M/E^/c queueing 
system can now be formally described. 

A. DESCRIPTION OF THEPROPOSED COMPUTATION METHOD 

After analysis of data and the decision has been made 
that the particular system under study can be approximated 
by M/E^/c system as described in Chapter II, and also having 
the estimators for X, k and 3 obtained (by maximum likeli- 
hood or method of moments) , one is ready for computations 
to get the desired information for the system. 

1. Delay Probability 

First compute service rate y as 




Then compute traffic intensity p and offered load a by 

X , X 

p = — and a = — 

cy y 

respectively. Then enter the chart in Appendix A with a to 
get delay probability P(T >0). Use the charts in [10] or 
[11] if number of servers c>16. 
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2 . Average Waiting Time and Average Queue Length 
First compute a system constant, namely b, by 



b 



2c y (1-p) 



1 



+ 



1 

k 



Then average waiting time W q 



is computed as follows 




P(T >0) 

^ mim 

b 



Also, average queue length L is given by 

'"‘d 




3 . Average Time in System and Average Number in System 
Average time in system W and average number in sys- 
tem L will be computed by 

W = W + - and L = L + - 

q y q y 



respectively. 

4 . Waiting Time Distribution 

The CDF of waiting times is given by 



F Tq (t) = l-P(T g >0)e bt , t>0 

where b is the system constant computed in step 2. 

Then the following probabilities can be computed 

(i) P(T q >t x ) = P(T q >0)e -bt l 

(ii) p ( t 1 <T q <t 2 ) = P (T q >0) (e _bt l - e~ bt 2) 

(iii) P(T <t.) = 1-P(T >0)e” bt l 

J. 

as needed by the user. 
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B. DISCUSSION 



1. Advantages of the Method 

(i) Needless to say, the most important advantage 
of the method is simplicity. The kind of information listed 
above can be computed in minutes given c, X, k and y which 
would be needed to compute anyway. All of the computations 
can be laid down on one regular size page so that it is 
very easy to follow by somebody else. 

(ii) In most cases, one has to compute the above 
listed information to determine the effect of changing the 
number of servers c on the selected system variables. This 
multiplies the savings of time and computational effort. 

(iii) No tables are necessary except for the delay 
probability chart. No interpolation would be needed. 

2 . Disadvantages of the Method 

(i) Since the method gives an approximation, it 

is not too precise even though the results are almost always 
in + 5 percent of the actual value, and in + 2 percent for 
most cases. 

(ii) The approximation of state probabilities with 
this method does not appear to be direct. 

3 . Applications 

One of the characteristic applications of the queue- 
ing theory is to investigate the effect of changing the 
parameters or number of servers on the measure of effective- 
ness (MOE) assigned by the management. This measure of 
effectiveness can be delay probability, average waiting time 
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or probability of waiting time exceeding time t etc. The 
marginal utility of c + 1st server will be the difference 
between the values of MOE's for c and c+1 servers. Then the 
decision criterion whether or not to hire the c+lst server 
is how his cost compares to his marginal utility. 

Another managerial objective may be to require, for 
example, the probability of waiting time exceeds 4 min., 

P (T >4) = 0.05. The number of servers necessary to achieve 
this purpose will then be of interest. 

It is possible to extend the variety of examples. 

The common property of the above examples is that precision 
greater than that of suggested approximation is not needed 
to make the decision. 

4 . Final Remarks 

The author wishes to emphasize that without the 
numerical tables provided by Hillier and Lo [6] , it would 
have been impossible to develop such an approximation method. 
The method may need some modifications as new tables for 
the cases not covered in [6] become available. 
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APPENDIX A 



CHARTS FOR DELAY PROBABILITY 

In the two following charts, the delay probabilities 
can be found for M/M/c system which were shown to be very 
close to those of M/E^/c system, with number of servers c 
up to about 16. One has to enter to chart with offered 
load a as abscissa, then delay probability P (T >0) is the 
ordinate value of the point on the proper c curve which 
has a as abscissa value. 



0 
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Offered Load 
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Offered Load 



APPENDIX B 



SAMPLE PAGES FROM NUMERICAL TABLES 



$utf 

I 



state 



9 I 
I I 
l I 
1 I 

4 I 

5 I 
9 I 
) I 

I I 



|TAfE 

l 



PIN- II 

0.B179210 

0. 16*1*78 

0.01666797 
0.0011828*2 
6.689l*4'-05 
3,236087' -06 
|.*00809'-07 
M • 1 • 



Pl*-I I 

0.6652663 
0.269*672 
0.05593797 
0.0082326*9 
0.0009823202 
0.000102*607 
0.799626'— 06 
8.890 805'-07 
7.853 90 8* -08 
M - 1 • ' 



P I N- t I 

0.3352*89 
0.329 5022 
0.1033992 
0.02*288*4 
0.00*626235 
0.0007867*96 
0.000125*372 
1.93*61 7 '-05 
2 . 9*0 7*5 ' -06 
M • t • 



STATE 

I 

0 I 

1 I 
I I 
3 I 

* I 

5 I 
9 I 



STATE 

I 

0 I 

1 I 

I I 
3 I 

* I 
5 I 

* I 

7 I 

8 I 

9 I 

10 I 

II I 



N - 



PIN<-I J 

0.8179210 
0.9820789 
0.9987*69 
0.9999297 
0.4999966 
0.9999 99 8 
0.94 99999 



P1N<-1 I 

0.6652663 
0.93*7536 
0.9906716 
0.99890*2 
0.9998868 
0.9999892 
0.9999990 
0.99 99999 
0.9999999 
• 8 



STATE 

1 

7 I 

8 I 

9 I 

10 I 

11 I 

12 I 

13 I 

C • 

STATE 

1 



PAN* I I 

5. 599922 * -09 
2.1 l 7**7' -Id 
7* 71 *0 W -1 2 
2.7*9*28 ' -13 
9.68) >72' -1 5 
5.595385* -16 
1. 190273**17 
2 • MO 



PI N<*1 I 

0.9999999 
0.9999999 
0.9999999 
0.4999999 
0 . 9999999 
0.4999999 
0.9999999 



P 4 N<- 1 1 



STATE 

( 

0 I 

1 I 

2 I 

3 I 



10 

11 



4 I 


6. SO 7781* -04 


0.4999994 


12 1 


10 1 


5.88745)'-10 


0.9994999 




11 1 


5.038657 • -l l 


0. 9949499 




12 I 


4. 4022l4'-12 


0. 9999999 


STATE 


13 1 


3.8122 38 • -l 3 


0.9999997 


1 


14 | 


5.30)589* -14 


1.000030 




13 | 


2.86)789*-l5 


1.000000 


0 


16 I 


2.482840' -16 


1.000330 


1 1 


17 1 


2. L 32 598 ' -1 7 


1.000000 


2 


C • 


2 . MO • 


0.30 


5 

4 



P 1 N<- 1 1 


1 


0.3352489 


9 


0.8647310 


10 


0.9701505 


11 


0.9944387 


12 


0.9990650 


13 


0.9998517 


14 


0.9999772 


15 


0.9994965 


16 


0.9999995 


17 



PIN-I 1 

O. *233869 
0.3532261 
0.1535*5* 
0.0 50 0 5 8 3 2 
0.01359216 
0.003360750 
0.0007926517 
0.0001831328 
*. 200506 *-05 
9.618*29* -06 
2. 202535* -06 
* - l • 



PIN-l ) 

0.3265101 

0.3*69793 

0.19728*6 

0.08362*06 

0.030*2698 

0.01025*93 

0.003339760 

0.00107373* 

0.0003*392*2 

0. 000 11009*0 
3.52*6*** -05 

1. !2857**-05 



PIN-1 ) 

1 0.2824307 

l 0.33*1386 

\ 0.21273*0 

| 0.1018312 

| 0.0*223033 

| 0.0163*158 

| 0.0061*1*69 

| 0.00228*530 

| 0.000 3*755*1 

| 0.00031*3163 

| 0.0001165765 

| *.J2*l22'-05 

M • 1 



STATE 

l 

0 I 

1 I 

2 I 

3 I 
* I 
6 I 

6 I 

7 I 

8 I 

9 I 

10 I 

11 I 

12 I 



PIN-l ) 

0.2*2259* 
.313*312 
.2231 736 
_. 1197532 
.036123*1 
j. 02*71 780 
0.01061 913 

0,00*52*238 
0.0014? J40* 
0.0003174**1 
0.0003*777*9 
0.0001*7377* 
6.2 3 3060 * —05 



0. 

0. 

0. 

0. 

0. 



> 8 * 


C - 2 


STATE 


P 1 N<- 1 1 


I 


0.4233869 


11 1 


0.7766131 


12 1 


0.9319585 


1 3 1 


0.9820168 


14 1 


0.99 56089 


15 1 


0.9989697 


16 1 


0.9997624 


17 | 


0, 9999*55 


IB 1 


0.99«9875 


19 I 


0.9949971 


20 1 


0.9999993 


21 1 


• 8 • 


C - l 


state 


P IN<- 1 > 


1 


0.3265101 


12 1 


0.6734899 


15 1 


0. 8707746 


14 | 


0.9543986 


15 1 


0.9848256 


16 1 


0.9950906 


17 I 


0. 998*204 


18 1 


0.9994941 


19 1 


0.9998381 


20 1 


0.9994*8 1 


21 1 


0.999983* 


22 1 


0.9 9 499 *6 


23 1 


• 8 . 


C • 
STATE 


pinou 


1 


0.2829307 


12 1 


0.6170692 


13 1 


0.8298033 


14 1 


0.9 316 5*6 


13 1 


0.9738831 


16 I 


0.9902267 


IT 1 


0.9963682 


18 1 


0.9986528 


19 1 


C. °995003 


20 1 


0. i*981*6 


21 1 


0.9999312 


22 I 


0.99997*3 


23 1 


• 8 


C • 

STATE 


P 1 N< — I 1 


1 


0.2*22594 


13 1 


0.557 7406 


14 | 


0.79091*2 


15 1 


0. 9 0066 7 5 


16 I 


0.93*7929 


17 1 


0.9815107 


18 1 


0.9921298 


19 


0.9966531 


20 


0. 99857 70 


21 


0 . *9039*'* 


22 


0.99 97*2 7 


23 


0.9998906 


2* 


0.999953* 


23 



PIN* 1 ) 

*. *48650 * -07 
6.724858'-08 
1.01 7040 ' -08 
1. 53 9943 * -09 
2.329 388* -10 
3.5262 66 * — ! 1 
5.3 3837 0 * -1 2 
8.08164** -1 3 
1.223*53 *-13 
AnO 



PIN-l) 

3.0*51*8* -07 
1.1558 76* -07 
2 .648*09 * -0 8 
6.063323* -09 
1, 3904*4 • -09 
).l 85936* -13 
7.299952* -11 
1.6726*3' ll 
3.832525* -12 
8.731*72* -13 
2.012101' -13 

» . KMO 



PIN-1 ) 

5.613901*-06 
1.15726J*-36 
3.7C5367'-07 
1.1*67 18* -07 
3.80018*' -08 
1.21 6919 • -08 
3.396843 * —39 
l.2*Te08* -09 
3.996063 '-13 
1.2 796*5 ' —l 0 
*.097760'-!* 
1. 31 22 10* -1 1 
2 * MO 



PIN* 1 I 

1.60 39 95 *-05 
3.9*9960' -36 
2. 2071 1 7* -06 
8.1»7226*-07 
5.0 t T 022 * -37 
1.126572*-0T 
*. 1789?«' -08 
1.5501 73' -08 
5.730316* -J9 
2,13 3059* -09 
7.912302* -10 
2.935tl2'-13 
2 1 MO 



P4NC-W 

0.99^9999 
0.9049999 
0.9999499 
0.9999949 
0*9999949 
0. 9999949 



7 

8 

9 

10 

11 

12 

13 



PIN-l ) 

0.20*2829 
0. 291*3*1 
0.2275953 
0.1 159851 
0.0713*6*9 
0.03559238 
0.01733*83 
0.003)877*8 
0.00*042620 
0.001937785 
0.0009*38477 
0.0004569 793 
0.0002207896 
4-1 • 



PlN?l) 

0. 1683107 
0.2623786 
0.22*9728 
0.1*88*35 
0.08739370 
0.0*579246 
0.02675223 
0.01*59134 
0.007949822 
0.004330918 
0.002359518 
0.001285322 
0.000 7003895 
0.0003815931 



PCN<* 1 1 

0.20*2879 
0.4957170 
0.723312* 
0.6592975 
0.9)08*91 
0.9664370 
0.96)7718 
0.9921 *95 
0.9962122 
0.9981 700 
0.4991158 
0.4995728 
0.9997936 
• 8 • 



STATE 

1 

13 

1* 

15 

16 

17 

18 

19 

20 
21 
22 
23 
2 * 
25 
C 



0.9999999 

0.4999999 




N - 1 • 


0.9499999 


STATE 




,40 


1 


PIN- 1 1 




0 1 


0.1336715 


PtN<M 1 


1 1 


0.2286372 




2 1 


0.2142 547 


0.9999948 


) 1 


0.15621*4 


0. 9949999 


4 1 


0.1018311 


0.9999999 


5 1 


0.065*4116 


0.4999999 


6 I 


0.03891755 


0.9999949 


7 1 


0.02377329 


0.9959994 


8 1 


0.01451033 


0. 9949499 


9 1 


0.008856043 


0.9999999 


10 1 


0.005*05270 


0.9999999 


U 1 


0.00)299163 


0.9999584 


12 1 


0.002013691 


0.9999999 


13 1 


| 0.001229035 


.50 


14 I 


0.0 00 7 5 0 1 903 


13 1 


| 0.000*578896 



P|N<-1 ) 

0. 1696107 
0. *311392 
0.6561621 
0.8050057 
0.892399* 
0. 9*t 1 919 

0.9679**1 

0.9825354 

0.990*€53 

0.99*8162 

0.9471757 

0.9984612 

0.9991616 

0.9995452 

• 8 t 



ptN<-n 

0. 1556713 
0. 364 32 36 
0.578583* 
0.7347978 
0.8366290 
0.9000702 
0.9389877 
0.5627610 
0.5772713 
0.9 8612 74 
0.9913327 
0.9948318 
0.9968435 
0 . 99 80 7 -8 
0.99882*6 
0.4992827 



l* 

15 

14 

17 

18 

19 

20 
21 
22 
23 
2 * 
23 



. NHO • 


0.63 


PIN-l A 


P IN<- l ) 


0.0001066 7*9 


0.9949002 


3.15*022' -06 


0.949951 8 


2. *901 76 • -03 


0. 9949767 


1.203135' -05 


0.9999687 


5.81296*' -06 


0.99T9445 


2.8083** * -06 


0.99*9973 


l. 3)6955* -06 


0.9999987 


6.5561*8 *-07 


0.99*9993 


3. 1 4 76 1 5 • -3 7 


0.9999997 


1. 5 304*0 * -0? 


0.9999998 


7 . 394345 ' -08 


0.9*99949 


3.572594*-0tt 


0.9999999 


l. 726106'-O8 


0.9999999 


2 • MU • 


J. 70 


PIN-l > 


P IN<- 1 1 


0.0002079035 


0.9947511 


0.0001132721 


0.9998644 


6.1 7 1*0* *-0 5 


0.94*9261 


3. 362367* -05 


0.4999597 


1.831919' -05 


0. 9999721 


9.980944 • —36 


0.9999880 


S-43?867'-06 


u, 9499934 


2.962714'— J6 


0.9999964 


l.6l4l73'-06 


0.4499980 


8.79451 5 ' -0 7 


0.9994489 


4. 791 517' -07 


0.49*9494 


2 .610562 * -3 7 


0.9949996 


1.422313* -07 


C. 9999198 


7. 7491 65* -08 


0.9449999 


2 , MO • 


0.75 



16 I 

17 I 

18 I 

19 I 

20 I 

21 I 

22 1 

23 I 
2* I 

25 I 

26 I 

27 I 

28 I 

24 I 



30 

31 



PlN-l > 

0.000279*797 

0.00017038*5 

0. 00010. 1187 

6. 3 550*1 * -05 
3.378892* -J5 
2. 367550 * -05 

1. **5362 ' —35 
8.8201*8* -J 6 

5.3 8 350? * -06 
3.285902* -Oo 
2.005598 * — 06 
1.22*l*6'-96 
7.471757'-07 
* .*60*96 • -0 7 

2. 78) 565 '-0 7 
1.698990* -07 



P1N2-I) 

0.4995621 
0.959732 7 
0.9998369 
O. 599900* 
0.9999)92 
0. 9999o2 9 
0.595977* 
0.9999862 
0.999941 5 
0.99V9”»*8 
0.9999468 
0. *449960 
0.9449483 
0.4459943 
0.9999993 
0.9949597 



UNO • 0 . so 



P 1 N<— I > 


0. 999 9983 


STATE 

1 


PIN-l 1 


P 1 NO I I 


3.8999794 


0.9999998 


0 1 


0.1047099 


0.10*7094 


0.9999999 


1 1 


0.1905802 


0. 2932901 


0.9999999 


2 I 


0.1943588 


0 • *8 >6*89 


0.9999999 


3 1 


0, 1553159 


0. 6*31628 


0.9999999 


4 1 


0.1120363 


0.7371992 


0.97 9999 9 


5 1 


0.07749480 


0.93*69*0 


0.9989999 


6 I 


0.03289866 


0.8875927 


0.9999999 


7 1 


0.03598627 


0.9255784 


0,5999999 


8 I 


0.02*46579 


0. 9 *80**7 


0.9999999 


9 1 


0.016632 78 


0.96*6775 


1. 55 


10 1 


0.01130789 


0.975935* 



PIN<-1 1 



9999905 
,9990945 
,9999987 
,9999995 
.9499998 
, 9959999 
,9999999 
.4994994 
0.9999999 
0. 9999999 
0.9999999 
0.9999949 



15 I 
1* I 



17 I 

18 I 



STATE 

1 

0 I 

1 I 

2 I 

3 I 



0 .00 7687 83 7 
0.003226709 
0.003353*71 
0.002*15890 
0 . 0016 * 2 * 8 * 
0.031116673 
). OCO 75 51 899 
). 0005161*87 
• 1 • 



0.96)6733 

0.9898999 

0.992*535 

0.99*869* 

0.9965119 

0.9976285 

0.9983677 

0.5989038 



STATE 

I 

19 

20 
21 
22 
23 
2* 

25 

26 
27 
23 

29 

30 

31 

32 I 

33 I 

34 | 

35 I 

36 | 

37 | 
C • 



2.67383)' 
1.136978’ 
4.83*706 
2.0539) *' 
6. 7*1893 
3. 7172*2 
1.3 836 6 9 
6.7<l)*2 
2.853 C 
1.213)29 
5.1678/0 
2.197502 
9.3*4)05 



-05 
-05 
-06 
-36 
-0 7 
-07 
-07 

> -j8 
1 -08 
1 -09 

> -09 
• -04 
*-lO 



PI N<" I I 

0.9999802 
0.9999913 
0.9999964 
0.9499995 
0.5999995 
0.9999597 
0.5999999 
0.9999999 
0.9449999 
0.9999944 
0.9994999 
0. 9999999 
fe. 9999999 



5 I 

6 I 

7 I 

8 I 

9 I 

10 I 

11 I 

12 I 

13 I 

14 I 

13 I 

16 I 

17 I 

18 I 

19 | 



21 

22 

23 

24 



PIN-tl 


PINO! 1 


STATE 

1 


0.07578546 


0.075785*6 


25 


0.1484290 


0.22*21*5 


26 


0.1641668 


0.338381* 


2? 


0. 1436229 


0. 532C042 


28 


0.1138887 


0.6*569)1 


29 


0.08706796 


0.7529611 


50 


0.06581450 


0. 7987756 


31 


0.0*961175 


0.8*8)87* 


32 


0.03738017 


0.8357673 


5) 


0.02316330 


0.9139310 


3* 


0.02121974 


0.9351507 


35 


0.01 593012 


0.9511389 


36 


0.3120463* 


0.96)183) 


37 


0.039076413 


0.9722617 


38 


0.0063 38637 


0.9791003 


34 


0.003152658 


0. 98*2330 


*0 


0.00 )392306 


0.9881333 


41 


0. C 0 2**2 51 51 


0.9910604 


42 


0.03720)976 


0.99)2644 


45 


0. GOI 06 O 6 OI 


0.99*9 250 


44 


0.001261 192 


0.9961762 


43 


0.000 >4271 92 


0.9971189 


*6 


0.000 7 10 298 3 


0.9973293 


47 


0.00)3351 79* 


0.998)6** 


48 


0.000-0323*5 


0.993 76 77 


49 



PIN-t ) 

0.0003509128 
0.0002385743 
0.0001621989 
0.0001102758 
7. 49 71 65 * -05 
3.097079* -05 
3. 465319* -05 
2. 355971 * HD5 
1.6017*8* -35 
l.008977*-05 
7. 403608* -06 
5.0 13 *74 '-06 
3.422096*-06 
2.326573*-06 
1.5 SI 761 *-06 
1.0 75388* -06 
7.3112 21* -0 7 
4.9 70663* -0 7 
3. 3 79593 *-C7 
2 . MO 



PIN-l I 

0^00010)8200 

0.000221*9153 

0.0001724778 

0. 0001299546 
9. 79l52?*-06 
7.377495'-V» 
5* 5 3 862 5 *-05 
*.lS8186'-03 

3. 133618* -05 
2,37 7623 *-0 3 
1.791437*-05 

1. 3 4 97 70 * -03 
1.31 6994 * -0 5 
7. 662516* -06 
5.77343l'-06 

4. 3 50048 * —06 
3.27 75 74'-0* 
2.469511 • -06 
1. B0O6 71 * —06 
1.40193J'-06 
1.056298 * —06 
7 • 9587 55 '-9 7 

5. 996530* -0 7 
4.3l8l66*-07 
3.4042V4' -07 



P1N2-1 ) 

,99923*0 
,999*9)3 
,9996535 
,999765* 
,9998*07 
,999891 7 
,9999244 
,9999499 
,9)99660 
,99997o8 
,949984 3 
0. 94V7P93 
0.4999927 
0.9999931 
0.9999966 
0.9999477 
0. 9999954 
0.9999989 
9999492 



0 

0* S3 



P I N<— 1 » 

,999071 3 
9993004 
,999*729 
,9996028 
,9997007 
,9997745 
W w 9998301 
0.9598720 
9944336 
9999273 
999945 2 
0.9999587 
0.9999684 
0.9994766 
0.4994*23 
0.9999866 
0.9959899 
0.9449924 
0. 9999943 
0.999995 7 
0.9 49946 7 
0.99*9*76 
0.9 999992 
0.9949986 
0.9449989 
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M * 1 


K - 8 . ♦ C 


= 2 




RHO 


P ( DELAY ) 


L ( GIVEN K) 


LQt GIVEN K) 


LQ FOR K~1 


RATIO 


0*10 


0.01792109 


0.2013269 


0.301326923 


0.002020202 


0.6568266 


0* 20 


0.06526637 


0.4105492 


0.01054924 


0.01666667 


0.6329545 


0*30 


3.1352489 


0.6365210 


0.03652108 


0.05934066 


0.6154479 


0.40 


0.2233869 


0.8917542 


0.09175420 


0.1523809 


0.6021371 


0.50 


0.3265101 


1.197244 


3.1972441 


0.3333333 


0.5917324 


0.55 


0.3829307 


1.380202 


0.2832032 


0*4770609 


0.5873531 


0.60 


0.4422594 


1.593808 


3.3938082 


0.6750000 


0.5834197 


0.65 


0.5042829 


1.851505 


0.551 5060 


0.9510822 


0.5798720 


o 

• 

o 


0.5688107 


2.175664 


0.7756644 


1.345098 


0.5766602 


0.75 


0.6356713 


2.60 6 503 


1.106503 


1.928571 


0.5737424 


o 

« 

CD 

o 


0. 7047099 


3.224415 


1.624415 


2.844444 


0.5710835 


0.85 


0.7757854 


4.216932 


2.516932 


4.426126 


0.5636536 


0.90 


0. 8487694 


6.146583 


4*3465 82 


7.673684 


0.5664272 


0.95 


0.9235438 


11.82589 


9.925892 


17.58717 


0.5643023 


0.98 


0.9692 222 


28.73332 


26.77332 


47.53494 


0.5632344 


0.99 


0.9845791 


56.86909 


54.88910 


97.51749 


0. 5628642 
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COMPUfER PROGRAM 



PREAMBLE 

PERMANENT ENTITIES 

THE SYSTEM OWNS A QUEUE 

EVERY SERVER HAS A STATE, AN IDENTI 

TEMPORARY ENTITIES 

EVERY JOB HAS AN ARR<» T I ME AND MAY BELONG TO THE QUEUE 
EVENT NOTICES INCLUDE ARRIVAL , END. OF . S I MULATI ON 
EVERY DEPARTURE HAS A CODE 

PRIORITY ORDER IS ARRI VAL , DEPARTURE , END. OF. SI MULATI ON 
NORMALLY, MODE IS REAL 

TALLY AV.Q.T AS THE MEAN AND VAR.Q AS THE VARIANCE OF R 
TALLY FREQ( 0 TO 15 BY 0.25) AS THE HISTOGRAM OF R 
TALLY AVER AS THE MEAN AND VARER AS THE VARIANCE OF ERLG 
ACCUMULATE MMM AS THE MEAN OF K 
ACCUMULATE KKK AS THE SUM OF MM 

ACCUMULATE FREStO TO 40 BY I) AS THE HISTOGRAM OF K 
DEFINE FREE TO MEAN 0 
DEFINE BUSY TO MEAN 1 

DEFINE X, Y, XX AND YY AS REA L, 1-DI MENS IONAL ARRAYS 

DEFINE F , Z AS REAL , 1-DI MENS I ONAL ARRAYS 

DEFINE NN,KK,STP,EPSTP,SToN,EP AS INTEGER VARIABLES 

DEFINE U, ERLG ,TOTi , AVG, MM AS REAL VARIABLES 

DEFINE STSl,SfS2,STS3 AS INTEGER VARIABLES 

DEFINE SD1,SD2,SD3 AS INTEGER VARIABLES 

DEFINE FLAG AS AN INTEGER VARIABLE 

DEFINE LAMBDA ,MU,S , TOTo QoTI ME , R » SUMSQR , D, BETA , ST , STT , START 0 T 
AS REAL VARIABLES 

DEFINE CUSTOMER, J,K, DEL AYER, AR,M,N, ALPHAS A2E.N0, I, EFLAG 
AS INTEGER VARIABLES 

END 



MAIN 

REAO CASE. NO 
READ J 

" J IS THE NUMBER OF SERVERS 
REAO ALPHA 

" ALPHA IS SHAPE PARAMETER OF THE SERVICE TIME DISTRIBUTION 
READ BETA 

" BETA IS SCALE PARAMETER OF THE SERVICE TIME DISTRIBUTION 
READ LAMBDA 

" LAMBDA IS ARRIVAL RATE 
READ EPSTP 

" EPSTP IS TOTAL NUMBER OF EPOCHS FOR THIS SIMULATION RUN 
READ FLAG 

" FLAG INDICATES WHICH PAIR OF ANTITHETIC VARIATES WILL BE 
' « USED 

READ SCI,SD2,SD3 
READ NN ,KK 

" NN IS DIMENSION OF ARRAY F 
" KK IS NUMBER OF ELEMENTS NE I IN ARRAY F 
RESERVE F ( * ) AS NN AND Z(*) AS NN 

" READ CDF TABLE OF THE CHI-SQR DISTN. WITH 2*ALPHA 
' ' DEGREES OF FREEDOM 

FOR 1=1 TO KK,DO 
READ F ( I ) 

LET F ( I )= 1. 0-F ( I ) * 

LOOP 

FOR 1=1 TO KK ,D0 
READ Z ( I > 

LOOP 

FOR I=KK+1 TO NN , DO 
LET F ( I )= 1.0 
LOO® 
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LET N.SERVER=J 
CREATE EACH SERVER 
LET K=0 
LET SIP =0 

LET STS1=SEED.V(SD1 ) 

LET STS2=SEED.V(SD2) 

LET STS3=SEED.V(SD3) 

START NEW PAGE 

SCHEDULE AN ARRIVAL NOW 

START SIMULATION 

STOP 

END 



EVENT ARRIVAL SAVING THE EVENT NOTICE 

DEFINE JJ AS AN INTEGER VARIABLE 

ADD 1 TO AR 

LET T=TIME.V*1440 

IF STP EQ 1 

DESTROY THIS ARRIVAL 

RETURN 

OTHERWISE 

ADD 1 TO K 

" GENERATE INTERARRIVAL TIME 
LET U=UNIFORM.F(OoO,loO,SDl ) 

IF FLAG EQ 1 
LET L=1„0-U 
ALWAYS 

LET AR«TIME=-LOG«E *F ( U ) /LAMBDA 
RESCHEDULE THIS ARRIVAL IN AR.TIME MINUTES 
LET JJ=0 

FOR EVERY SERVER, DO 

IF STATE( SERVER) =BUSY 

ADD 1 TO JJ 

ALWAYS 

LOOP 

IF JJ=0 

LET E FLAG=1 
CALL EPOCH 
GO TO AA 
ALWAYS 
LET J1=J-1 
IF JJ EQ J1 
LET MM=1.0 
ALWAYS 
• AA • 

•NEXT' FGR EVERY SERVER, DO 
IF STATE ( SERVER ) =FREE 
LET STATE (SERVER )=BUSY 
LET R=0<>0 

•ADO' ADD 1 TO CUSTOMER 
' • GENERATE SERVICE TIME 
LET U=UNI FORM«F(OoO,loO,SD2 ) 

IF FLAG EQ 1 

LET U=1.0-U 

ALWAYS 

CALL ERLNG 

LET DEP.TIME=ERLG 

LET IDENTI(SERVER)=DEPoTIM5 

SCHEDULE A DEPARTURE GIVEN DEPoTIME IN DEPoTIME MINUTES 

RETURN 

OTHERWISE 

LOOP 

CREATE A JOB 

LET ARRoTIME< JOB ) =T IME » V*1440 
FILE THIS JOB IN THE QUEUE 
RETURN 
END 
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EVENT DEPARTURE SAVING THE EVENT NOTICE 
DEFINE JJ AS AN INTEGER VARIABLE 
LET T=TIME.V*1440 
IF STP EQ 1 

DESTROY THIS DEPARTURE 

RETURN 

OTHERWISE 

SUBTRACT 1 FROM K 

LET JJ=0 

FOR EACH SERVER, DO 

IF STATE(SERVER)=BUSY 
ADD 1 TO JJ 
ALWAYS 
LOOP 

IF JJ EC 1 
LET EFLAG=0 
CALL EPOCH 
GO TO AA 
OTHERWISE 
IF JJ EQ J 

I F N* QUEUE EQ 0 
LET MM=0® 0 
ALWAYS 
ALWAYS 



• AA ' 

FOR EVERY SERVER, DO 

I F I DENT I ( SERVER ) =CODE (DEPARTURE ) 

IF THE QUEUE IS EMPTY 

LET STATE ( SERVER) 3 FREE 
DESTROY THE DEPARTURE 
RETURN 
OTHERWISE 

LET A,T=ARR.TIME(FoQUEUE ) 

REMOVE FIRST JOB FROM THE QUEUE 

DESTROY THIS JOB 

LET R=T-AoT 

ADD 1 TO CUSTOMER 

ADD I TC DELAYER 

ADD R TO TOTo Q«TIME 

ADD R**2 TO SUMSQR 

• • GENERATE SERVICE TIME 

LET U=UNI FORM •F(0o0,l«0,SD3) 

IF FLAG EQ I 

LET 1= 1 oO~U 

ALWAYS 

CALL ERLNG 

LET DEPoTIME=ERLG 

LET IDENTKSERVER) =DEPoTIME 

SCHEDULE A DEPARTURE GIVEN DEP.TIME IN DEP.TIME MINUSES 

RETURN 

OTHERWISE 

LOOP 

END 



EVENT END.CFo SIMULATION 

DEFINE FR, LUCKY AS INTEGER VARIABLES 

DEFINE DELI AS AN INTEGER VARIABLE 

RESERVE X(*) AS 60, Y{*) AS 60 AND YY(*) AS 60 

RESERVE XX(*) AS 60 

LET S=TIMEoV*1440 

START NEW PAGE 

LET MU=1.0/(ALPHA*BETA ) 

" MU IS SERVICE RATE 
LET G=L AMBDA/ MU 
• • G IS OFFERED LOAD 
LET RO=G/J 
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" RO IS TRAFFIC INTENSITY 

" COMPUTE AVG.WAIT .TIME FOR EQUIVALENT M/M/C QUEUE 

LET JFAC=I 

FOR 1=1 TO J-1 t DO 

LET JFAC=JFAC*(I+1) 

LOOP 

LET C=( <J*RO) **J)/(JFAC*< loO-RO) ) 

LET B=1.0 

FOR K = 1 TO J-l» DO 

LET KFAC=1 

FOR L=1 TO K-l , DO 
LET KFAC=KFAC*(L+1) 

LOOP 

ADD ( ( J*RO)**K)/KFAC TO B 
LOOP 

LET PON£=B+C 

“ PN0T=1.0/ (B+C) 

E = PNOT *( <J*RO)**J)/(J*JFAC*< ( loO-RO ) *#2 ) *MU ) 
VSQR=loO/ ALPHA 



LET 

LET 

LET 

LET 

LET 

LET 



V=SQRT.F< VSQR) 

W= ( E/2o0 ) # ( 1 oO+VSQR ) 
PMMC=E*J*MU* (loO-RO ) 



PRINT 1 LINE WITH CASE. NO 
CASE NO=**** 
SKIP 2 LINES 

PRINT 2 LINES AS FOLLOWS 
INPUT INFORMATION: 



AND FLAG AS FOLLOWS 
FLAG=* 



5RTP 2 CTFiK 

PRINT 5 LINES WITH STS 1 , SEED. V< SD1 ), STS2 , SEED. V <S D2 ) 
STS3,ScED.V(SD3) AS FOLLOWS 

STARTING SEEDS: FINAL SEEDS: 

j ”5555555557 5555555555 



II. ####*•**#$* 

III. ###****##* sjc #$*#*$){£ 

SKIP 2 LINES 

PRINT 1 LINE WITH ALPHA AMD BETA AS FOLLOWS 
ALPHAs#** BET A=**.**** 

SKIP 2 LINES 

PRINT 1 LINE WITH LAMBDA AS FOLLOWS 
L AMBDA=*##o**** 

SKIP 2 LINES 

PRINT 1 LINE WITH MU AS FOLLOWS 

MU=#**o#*** 

SKIP 2 LINES 

PRINT 1 LINE WITH RO AS FOLLOWS 
TRAFFIC INTENSITY, RO=*.**** 

SKIP 2 LINES 

PRINT 1 LINE WITH J AS FOLLOWS 
NUMBER OF SERVERS=** 

SKIP 2 LINES 

PRINT 1 LINE WITH G AS FOLLOWS 
OFFERED LOAD=**.**** 

SKIP 2 LINES 

PRINT 1 LINE WITH V AS FOLLOWS 
COEFFICIENT OF VARIATIONS**.*#** 

START NEW PAGE 
IF EFLAG=1 

ADD S-START.T TO SUMP 
ALWAYS 

LET LUCKY=CUSTOM5R-DELAYER 
LET PD= DELAYER/CUSTOMER 
PRINT 2 LINES AS FOLLOWS 
OUTPUT INFORMATION: 



» 



5kTI3-2-lTFJE3 

PRINT 1 LINE WITH AV.Q.T AS FOLLOWS 

MEAN QUEUEING TIME FROM S I MULA T ION = ***o #*** 

SKIP 2 LINES 

PRINT 1 LINE WITH W AS FOLLOWS 

COMPUTED MEAN QUEUEING TIME=***.**** MINUTES 
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, Qo T AS FOLLOWS 



SKIP 2 LINES 
PRINT i LINE WITH W-AV, 

DIFFERENCE®**.**** 

SKIP 2 LINES 
PRINT 1 LINE WITH VAR.Q AS FOLLOWS 
VARIANCE OF QUEUEING-TI ME=***o**** 

SKIP 2 LINES 

PRINT 1 LINE WITH MMM AS FOLLOWS 
AVERAGE QUEUE LENGTH®**.**** 

SKIP 2 LINES 

_ WITH CUSTOMER AS FOLLOWS 

SERVED IN SIMULATION PERIOD®***** 



WITH DELAYER 
OF CUSTOMERS 



AS FOLLOWS 
WITH DELAY 



ivis 



AS FOLLOWS 



FOLLOWS 

MINUTES 



PRINT 1 LINE 
TOTAL NUMBER 
SKIP 2 LINES 
PRINT 1 LINE 
TOTAL NUMBER 
LET PER IOD=S-STT 
SKIP 2 LINES 

PRINT i LINE WITH PD AS FOLLOWS 
PRO 8(W>0)=*.* ***** 

SKIP 2 LINES 
PRINT 1 LINE WITH PMMC 
PMMC=*. ****** 

SKIP 2 LINES 

PRINT 1 LINE WITH PERIOD AS 
SIMULATION TI ME=******. **** 

SKIP 2 LINES 

PRINT I LINE WITH KKK*1440 AS FOLLOWS 
TOTAL BUSY PERIOD=********o **** MINUTES 
FOR 1=1 TO 60,00 
LET XX( I ) =FREQ( I ) 

LOOP 

SUBTRACT LUCKY FROM XX(1) 

" COMPUTE AND PRINT F SUB I 

FOR 1=1 TO N , DO 

ADD XX ( I > TO FR 

LET Y(I )=1.0- ( FR/DELAYER ) 

IF FLAG EQ 0 
IF FR EC DELAYER 
LET N = I - 1 
GO TO REG 
OTHERWISE 
ALWAYS 

LET Y ( I ) = LOG.E.F(Y(I ) ) 

LET X(I )= 1/4. 0 

LOOP 

•REG' 

LET CUMFREQ=0 
START NEW PAGE 

" COMPUTE AND PRINT STATE PROBABILITIES 
PRINT 2 LINES AS FOLLOWS 
I X ( I ) FREQ ( I ) 



GREATER THAN ZERO=***** 



K5R 1=1 TO N,D3 
ADD XX( I > TO CUMFREQ 
SKIP 1 LINE 

PRINT 1 LINE WITH I,X(I),XX(I) AND CUMFREQ 
** **.*** ********* 

LOOP 

IF FLAG EQ 1 

CALL REGRESSION 
Al WAYR 

START NEW PAGE 
PRINT 2 LINES AS FOLLOWS 
N T ( N ) P (N ) 



AS 



CUMFREQ 



FOLLOWS 

********* 



FOE 1=1 

PRINT 1 

** 

LOOP 

RETURN 

END 



TO 40733" 
LINE WITH 



I,FRES(I)*1440o0 AND < FRES ( I ) * 1 440. 
AS FOLLOWS 

*****.**** * # ******** 



) /PERI 0 
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ROUTINE EPOCH 

DEFINE FI.N AS AN INTEGER VARIABLE 

LET T=TIME.V*l440 

IF EFLAG EQUALS I 

•• EPOCH STARTS 

LET START. T=T 

LET ST oN'= CUSTOMER 

LET TOTl=TOTo Q.TIME 

ADD 1 TO EP 

RETURN 

OTHERWISE v 

• ' EPOCH ENDS 
LET FINISH.T = T 
LET FI.N=CUSTOMER 
LET T0T2=T0T. Q.TIME 

LET AVG=( TQT2-T071 ) / ( ( F I. N-ST.N ) *1.0 ) 

LET PERICD=FINISHoT-STARToT 
LET CUM.AVG=T0T2/FI.N 
SKI 0 1 LINE 

LIST EP » START ..T, FINISH. T, ST. N,F I. N, TOT 1, TOT 2 
LIST PERI0D,T0T2-T0Tl,FIoN-SToN,AVG,CUM.AVG 
" CHECK FOR END CF SIMULATION 
IF EP EQ EPSTP 
LET ST P=1 

SCHEDULE AN END.OFo S I MULATI ON NOW 
ALWAYS 
RETURN 
END 



ROUTINE ERLNG 

DEFINE RL »LH, RH AS INTEGER VARIABLES 
IF U LT F ( 1 ) 

LET CHSQR=Z(I)*(U/F(1) ) 

GO TO ERLANG 
ALWAYS 

•• START BISECTION SEARCH 
IF U GT F (KK ) 

LET CHSQR=Z(KK)+<Z(KK)-Z{KK-l))*<U-F(KK) ) 

GO TO ERLANG 
ALWAYS 
LET LH=Q 
LET RH=NN 

•AA' LET RL = ( LH+RH > /2 
IF U LT F { RL ) 

LET RH=RL 
GO TO BB 
ALWAYS 
LET LH=RL 
•BB' IF RH=LH+1 

GO TO QUIT 
ALWAYS 
GO TO AA 
'QUIT* 

" DO LINEAR INTERPOLATION 

LET CHSQR=Z ( LH ) + ( ( U-F (LH) )*(Z(RH)-Z(LH) ) )/ ( F { RH ) -F( LH ) ) 
•ERLANG* 

LET ERLG=(BETA*CHSQR)/2.0 

RETURN 

END 
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